---
title: "Replication Archive: Ethnicity and Federalism Preferences in Ethiopia"
author: "Davis & Tellez"
date: today
format:
  html:
    toc: true
    toc-depth: 3
    code-fold: true
    embed-resources: true
execute:
  warning: false
  message: false
---

```{r setup}
#| label: setup

# Core libraries
library(tidyverse)
library(stargazer)
library(scales)
library(cregg)
library(broom)
library(ggrepel)
library(patchwork)
library(haven)
library(ggbump)
library(ggeffects)
library(pwr)
library(MetBrewer)

# Custom theme (from R/00-funcs.R)
theme_nice <- function() {
  theme_minimal(base_size = 14) +
    theme(panel.grid.minor = element_blank(),
          plot.background = element_rect(fill = "white", color = NA),
          plot.title = element_text(face = "bold"),
          axis.title = element_text(face = "bold"),
          strip.text = element_text(face = "bold", size = rel(0.8), hjust = 0),
          strip.background = element_rect(fill = "grey80", color = NA),
          legend.title = element_text(face = "bold"))
}

theme_set(theme_nice())

# Palette
palette <- MetBrewer::met.brewer(name = "Cross")
blue <- palette[8]
green <- palette[5]
red <- palette[2]

# Export helper (saves to replication/figures/)
ggsave_both <- function(x = last_plot(), folder = "figures",
                        name = "plot", type = "both",
                        width = 7, height = 7) {
  dir.create(folder, showWarnings = FALSE, recursive = TRUE)
  if (type %in% c("pdf", "both")) {
    ggsave(plot = x, filename = paste0(folder, "/", name, ".pdf"),
           width = width, height = height, dpi = 700,
           device = cairo_pdf)
  }
  if (type %in% c("png", "both")) {
    ggsave(plot = x, filename = paste0(folder, "/", name, ".png"),
           width = width, height = height, dpi = 700,
           device = "png")
  }
}

# Conjoint feature/level labels (used by script 06)
features <- tibble::tribble(
  ~feature, ~clean_feature,
  "hiring", "Who decides public hiring?",
  "police", "One national police or state police forces?",
  "culture", "Who decides official language, teaching of history?",
  "borders", "How are state borders drawn?"
)

levels <- tibble::tribble(
  ~level, ~clean_level,
  "Central government hires all civil servants", "The central gov't",
  "Each state hires own civil servants", "The states",
  "Each state has own police force", "State police",
  "One national police force", "National police",
  "Each state decides official language, teaching of history", "The states",
  "The central government decides official language, teaching of history", "The central gov't",
  "State borders are drawn according to geography", "Geography",
  "State borders are drawn according to majority ethnic group in the area", "Ethnicity"
)
```

# Main Text

## Figure 1: Ethnic Power Relations & Woreda Demographics

Source: `09_EPR figures.R`. Combines a bump chart of ethnic power relations over time (EPR Dataset) with woreda-level demographic distributions (2007 census).

```{r fig1-woreda-bump-combined}
#| label: fig-woreda-bump-combined
#| fig-width: 12
#| fig-height: 12

# EPR data
epr <- read.csv("data/EPR-2023.csv")
ethiopia <- epr |>
  filter(statename == "Ethiopia",
         group %in% c("Amhara", "Oroma", "Tigry",
                       "Other Southern Nations", "Somali (Ogaden)"))

# Census data
regions <- readxl::read_excel("data/ethphcensusatlas_data.xls", sheet = 2)
groups <- readxl::read_excel("data/ethphcensusatlas_data.xls", sheet = 7)
merged <- left_join(regions, groups)

wide <- merged |>
  select(R_NAME, W_NAME, ORO_PRCT, AMH_PRCT, TIG_PRCT) |>
  drop_na() |>
  pivot_longer(-c(R_NAME, W_NAME)) |>
  mutate(hilite = case_when(R_NAME == "Tigray" & name == "TIG_PRCT" ~ TRUE,
                            R_NAME == "Amhara" & name == "AMH_PRCT" ~ TRUE,
                            R_NAME == "Oromiya" & name == "ORO_PRCT" ~ TRUE,
                            TRUE ~ FALSE)) |>
  mutate(name = case_when(name == "AMH_PRCT" ~ "Amhara",
                          name == "ORO_PRCT" ~ "Oromo",
                          name == "TIG_PRCT" ~ "Tigray"))

# Bump plot
ethiopia <- ethiopia %>%
  group_by(from) %>%
  mutate(status = fct_rev(factor(status, levels = c("DOMINANT", "SENIOR PARTNER",
                                                     "JUNIOR PARTNER", "POWERLESS",
                                                     "DISCRIMINATED"))),
         status_numeric = as.numeric(status),
         status_offset = case_when(group == "Amhara" ~ status_numeric - .06,
                                   group == "Oroma" ~ status_numeric,
                                   group == "Other Southern Nations" ~ status_numeric - .03,
                                   group == "Somali (Ogaden)" ~ status_numeric + .06,
                                   group == "Tigry" ~ status_numeric - .06)) |>
  mutate(period = fct_inorder(paste0(from, "-\n", to)))

p1 <- ggplot(ethiopia, aes(x = period, y = status_offset, color = group, group = group)) +
  geom_point(size = 4) +
  geom_bump(size = 2, smooth = 8) +
  theme(axis.ticks = element_blank(),
        legend.position = "top") +
  scale_y_continuous(
    name = NULL,
    breaks = 1:5,
    labels = levels(ethiopia$status),
    sec.axis = sec_axis(~ .,
                        labels = levels(ethiopia$status),
                        name = NULL)) +
  labs(x = NULL, y = NULL, color = NULL, shape = NULL,
       title = "Changing ethnic power relations in Ethiopia",
       subtitle = "Source: Ethnic Power Relations Dataset, 2023.") +
  theme(panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank()) +
  scale_color_manual(values = palette[c(2, 4, 6, 8, 10)])

p2 <- ggplot(wide, aes(y = fct_rev(R_NAME), x = value, fill = hilite)) +
  ggbeeswarm::geom_quasirandom(shape = 21, alpha = 0.5,
                                varwidth = FALSE,
                                bandwidth = 0.9) +
  facet_wrap(vars(name)) +
  scale_fill_manual(values = c("gray20", red)) +
  theme(legend.position = "none") +
  labs(x = "Percent of residents in woreda belonging to ethnic group", y = "Region",
       title = "Ethiopia's demographic composition",
       subtitle = "Source: 2007 census. Each point represents one woreda.") +
  scale_x_continuous(labels = scales::percent_format(scale = 1)) +
  scale_y_discrete(labels = scales::label_wrap(15))

p1 / p2 + plot_annotation(tag_levels = "A")

ggsave_both(name = "woreda-bump-combined", width = 12, height = 12)
```

## Figure 2: Federalism Vignette Effects

Source: `05-vignette-experiment.R`. Treatment effects of ethnic composition information on support for federal power.

```{r fig2-federalism-vignette}
#| label: fig-federalism-vignette
#| fig-width: 7
#| fig-height: 7

baseline <- read_rds("data/clean_baseline.rds") |>
  mutate(federalism_treat = fct_relevel(federalism_treat,
                                        "Control (no info)",
                                        "R is in majority",
                                        "R is in minority",
                                        "R is in majority, but most co-ethnics in minority")) |>
  mutate(is_oromo = ifelse(q9 == "Oromo", 1, 0))

mod <- lm(q45_1 ~ federalism_treat, data = baseline)

preds <- ggeffects::ggeffect(model = mod, terms = "federalism_treat") %>%
  as_tibble()

ggplot(preds, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(fatten = 4, linewidth = 3, shape = 21, fill = "white", color = blue) +
  scale_y_continuous(limits = c(4, 9)) +
  scale_x_discrete(labels = scales::label_wrap(20)) +
  labs(x = "Experimental conditions: in newly defined state...",
       y = "All political power is given to... \n(0 = states, 10 = central gov't)")

ggsave_both(name = "federalism_vignette_effects")
```

## Figure 3: Heterogeneous Effects (One-Party & Tolerance)

Source: `05-vignette-experiment.R`. Combined interaction plots for one-party support and ethnic tolerance moderators.

```{r fig3-federalism-interact}
#| label: fig-federalism-interact
#| fig-width: 12
#| fig-height: 10

# One-party rule interaction
mod_1p <- lm(q45_1 ~ federalism_treat * one_party, data = baseline)

preds_1p <- ggeffects::ggpredict(mod_1p, terms = c("federalism_treat", "one_party [1, 3, 5]"))

party_plot <- ggplot(preds_1p, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, color = group)) +
  geom_pointrange(fatten = 4, linewidth = 3, shape = 21, fill = "white",
                  position = position_dodge(width = 0.3)) +
  scale_y_continuous(limits = c(4, 9)) +
  scale_x_discrete(labels = scales::label_wrap(20)) +
  labs(x = "Experimental conditions: in newly defined state...",
       y = "All political power is given to... \n(0 = states, 10 = central gov't)",
       color = "Only one party allowed to run for elections:") +
  scale_color_manual(labels = c("1" = "Strongly disagree", "3" = "Neither agree nor disagree", "5" = "Strongly agree"),
                     values = palette[6:8]) +
  theme(legend.position = "top")

# Tolerance interaction
mod_tol <- lm(q45_1 ~ federalism_treat * tolerance, data = baseline)

preds_tol <- ggeffects::ggpredict(mod_tol, terms = c("federalism_treat", "tolerance [1, 3, 5]"))

tol_plot <- ggplot(preds_tol, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, color = group)) +
  geom_pointrange(fatten = 4, linewidth = 3, shape = 21, fill = "white",
                  position = position_dodge(width = 0.3)) +
  scale_y_continuous(limits = c(4, 9)) +
  scale_x_discrete(labels = scales::label_wrap(20)) +
  labs(x = "Experimental conditions: in newly defined state...",
       y = "All political power is given to... \n(0 = states, 10 = central gov't)",
       color = "Feelings towards \noutgroup neighbors") +
  scale_color_manual(labels = c("1" = "Stong dislike", "3" = "Neither like or dislike", "5" = "Strong like"),
                     values = palette[6:8]) +
  theme(legend.position = "top")

party_plot / tol_plot + plot_annotation(tag_levels = "A",
                                        tag_prefix = "(",
                                        tag_suffix = ")")

ggsave_both(name = "federalism_interact_combine", width = 12, height = 10)
```

## Figure 4: Conjoint Marginal Means

Source: `06-decentralize-cjt.R`. Marginal mean estimates from the decentralization conjoint experiment.

```{r fig4-conjoint-mms}
#| label: fig-conjoint-mms
#| fig-width: 7
#| fig-height: 7

fed_cjt <- read_rds("data/clean-decent-conjoint.rds")

form <- outcome ~ hiring + police + culture + borders

mms <- cj(data = fed_cjt,
           formula = form,
           id = ~ r_id,
           estimate = "mm") %>%
  mutate(
    se = (upper - lower) / (2 * 1.96),
    z_value = estimate / se,
    p_value = 2 * (1 - pnorm(abs(z_value)))
  )

mms %>%
  left_join(features) %>%
  left_join(levels) %>%
  ggplot(
    aes(x = estimate, y = clean_level,
        xmin = lower, xmax = upper,
        color = clean_feature)) +
  geom_pointrange(fatten = 3, linewidth = 2, shape = 21, fill = "white") +
  theme(legend.position = "none", strip.text = element_text(hjust = 0.5)) +
  labs(x = "Marginal mean estimates\nPr(Profile chosen | attribute)",
       y = NULL) +
  facet_wrap(vars(clean_feature), ncol = 1, scales = "free_y") +
  scale_x_continuous(labels = scales::percent) +
  scale_y_discrete(labels = scales::label_wrap(30)) +
  scale_color_manual(values = palette[c(2, 4, 6, 8)]) +
  geom_vline(xintercept = .50, lty = 2)

ggsave_both(name = "federal_cjt_mms")
```

## Figure 5: Conjoint by Ethnic Violence Support

Source: `06-decentralize-cjt.R`. Marginal means split by support for ethnic violence.

```{r fig5-conjoint-eth-violence}
#| label: fig-conjoint-eth-violence
#| fig-width: 8
#| fig-height: 7

fed_cjt <- fed_cjt %>%
  mutate(ev_dichotomous = factor(ev_dichotomous),
         three_ethn = factor(three_ethn))

mm_eth_violence <- fed_cjt %>%
  select(outcome, hiring, police, culture, borders, ev_dichotomous) %>%
  drop_na() %>%
  cj(data = .,
     formula = outcome ~ hiring + police + culture + borders,
     by = ~ ev_dichotomous,
     estimate = "mm") %>%
  mutate(ev_dichotomous = factor(ev_dichotomous, levels = c(0, 1),
                                 labels = c("No", "Yes")))

mm_eth_violence |>
  left_join(features) %>%
  left_join(levels) %>%
  ggplot(aes(x = estimate, y = clean_level,
             xmin = lower, xmax = upper,
             color = ev_dichotomous)) +
  geom_pointrange(fatten = 3, linewidth = 2, shape = 21, fill = "white",
                  position = position_dodge(width = 0.5)) +
  theme(legend.position = "top", strip.text = element_text(hjust = 0.5)) +
  labs(x = "Marginal Means\nPr(Profile chosen | attribute)",
       y = NULL,
       color = "Any Support for Ethnic Violence?") +
  facet_wrap(vars(clean_feature), ncol = 1, scales = "free_y") +
  scale_x_continuous(labels = scales::percent) +
  scale_y_discrete(labels = scales::label_wrap(20)) +
  scale_color_manual(values = c(blue, red)) +
  geom_vline(xintercept = 0.5, lty = 2)

ggsave_both(name = "fed_cjt_eth_violence", width = 8)
```

## Table 1: Conjoint Attributes

Table 1 is a static description of the conjoint design attributes and levels. It is not generated programmatically.

| Attribute | Level 1 | Level 2 |
|-----------|---------|---------|
| Who decides public hiring? | Central government hires all civil servants | Each state hires own civil servants |
| One national police or state police forces? | One national police force | Each state has own police force |
| Who decides official language, teaching of history? | The central government decides | Each state decides |
| How are state borders drawn? | According to geography | According to majority ethnic group |

---

# Appendix

## The Population and Sample

### Figure A1: Government Employment by Education

Source: `07_afrobarometer group comparison.R`. Government employment rates across education levels from Afrobarometer.

```{r figA1-govt-education}
#| label: fig-govt-education
#| fig-width: 10
#| fig-height: 8

ab <- haven::read_sav("data/afrobarometer.sav")

clean_ab <- ab |>
  select(Q97, Q95D, Q11A, Q11B, Q11C, withinwt_hh, Q13, Q55A, Q55B, Q55C, Q55D, Q55E, Q1) |>
  mutate(Q97 = as_factor(Q97),
         Q95D = as_factor(Q95D),
         Q13 = as_factor(Q13)) |>
  mutate(govt = case_when(Q95D == "Government" ~ 1,
                          Q95D == "Missing" ~ NA_integer_,
                          Q95D == "Refused" ~ NA_integer_,
                          Q95D == "Don't know" ~ NA_integer_,
                          TRUE ~ 0)) |>
  mutate(college = case_when(Q97 %in% c("Some university", "University completed", "Post-graduate") ~ 1,
                             Q97 %in% c("Don't know", "Refused", "Missing") ~ NA_integer_,
                             TRUE ~ 0)) |>
  mutate(vote = case_when(Q13 == "I voted in the election" ~ 1,
                          Q13 %in% c("Missing", "Refused", "Don't know") ~ NA_integer_,
                          TRUE ~ 0)) |>
  mutate(across(c(Q11A, Q11B, Q11C), ~ ifelse(. > 4, NA, .))) %>%
  mutate(across(c(Q11A, Q11B, Q11C), ~ ifelse(. == -1, NA, .))) %>%
  mutate(participation = Q11A + Q11B + Q11C) %>%
  mutate(protest = case_when(Q11C == 0 ~ 0,
                             Q11C == 1 ~ 0,
                             Q11C == 2 ~ 1,
                             Q11C == 3 ~ 1,
                             Q11C == 4 ~ 1)) %>%
  mutate(across(c(Q55A, Q55B, Q55C, Q55D, Q55E), ~ ifelse(. > 4, NA, .))) %>%
  mutate(across(c(Q55A, Q55B, Q55C, Q55D, Q55E), ~ ifelse(. == -1, NA, .))) %>%
  mutate(media_eng = Q55A + Q55B + Q55C + Q55D + Q55E) %>%
  mutate(young_25 = case_when(Q1 > 25 ~ 0,
                              Q1 <= 25 ~ 1))

clean_ab |>
  group_by(Q97) |>
  summarise(govt = weighted.mean(govt, w = withinwt_hh, na.rm = TRUE)) |>
  filter(Q97 != "Don't know") |>
  ggplot(aes(y = Q97, x = govt)) +
  geom_col() +
  scale_y_discrete(labels = scales::label_wrap(30)) +
  labs(title = "Government employment by education level",
       y = "Highest education level completed",
       x = "% Employed by the government",
       subtitle = "Source: Afrobarometer, Ethiopia (2020).",
       caption = "Note: employment rates use sampling weights.") +
  theme_nice() +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1))

ggsave_both(name = "govt_education", width = 10, height = 8)
```

### Figure A2: Participation Comparison (Afrobarometer vs. Study)

Source: `08_compare_afro_baseline.R`. Compares protest and community meeting attendance between Afrobarometer and study sample.

```{r figA2-participation-comparison}
#| label: fig-participation-comparison
#| fig-width: 10
#| fig-height: 6

baseline_comp <- read_rds("data/clean_baseline.rds") |>
  janitor::clean_names()

afro_comp <- haven::read_sav("data/afrobarometer.sav") |>
  janitor::clean_names()

# Recode Afrobarometer
afro_comp <- afro_comp |>
  mutate(across(c(q11a, q11c), ~ ifelse(. > 4, NA, .))) |>
  mutate(across(c(q11a, q11c), ~ ifelse(. == -1, NA, .))) |>
  mutate(q11c = case_when(q11c == 0 ~ 0, q11c == 1 ~ 0,
                          q11c == 2 ~ 1, q11c == 3 ~ 1, q11c == 4 ~ 1)) |>
  mutate(q11a = case_when(q11a == 0 ~ 0, q11a == 1 ~ 0,
                          q11a == 2 ~ 1, q11a == 3 ~ 1, q11a == 4 ~ 1)) |>
  select(withinwt_hh, radio = q55a, tv = q55b, newspaper = q55c,
         attend_comm = q11a, attend_protest = q11c,
         one_party = q20a, compromise = q31, diverse = q68,
         trust = q83, federalism = q76c_eth, unfairly_treated = q84c) |>
  mutate(across(-withinwt_hh, haven::as_factor)) |>
  mutate(attend_comm = as.double(attend_comm)) |>
  mutate(attend_comm = case_when(attend_comm == 1 ~ 0, attend_comm == 2 ~ 1,
                                 TRUE ~ as.double(NA))) |>
  mutate(attend_protest = as.double(attend_protest)) |>
  mutate(attend_protest = case_when(attend_protest == 1 ~ 0, attend_protest == 2 ~ 1,
                                    TRUE ~ as.double(NA)))

baseline_comp <- baseline_comp |>
  mutate(q13_2 = case_when(q13_2 == "Never" ~ 0,
                           q13_2 == "Once or Twice" ~ 1,
                           q13_2 == "More than twice" ~ 1,
                           q13_2 == "More than 5 times" ~ 1,
                           q13_2 == "More than 10 times" ~ 1)) |>
  mutate(q13_1 = case_when(q13_1 == "Never" ~ 0,
                           q13_1 == "Once or Twice" ~ 1,
                           q13_1 == "More than twice" ~ 1,
                           q13_1 == "More than 5 times" ~ 1,
                           q13_1 == "More than 10 times" ~ 1)) |>
  select(radio = q11_1, tv = q11_2, newspaper = q11_3,
         attend_comm = q13_1, attend_protest = q13_2,
         one_party = q18, compromise = q19, diverse = q20,
         trust = q21, federalism = q23, unfairly_treated = q24) |>
  mutate(compromise = factor(compromise, levels = c("I agree strongly with A",
                                                    "I agree with A",
                                                    "I agree with B",
                                                    "I agree strongly with B",
                                                    "I agree with neither"))) |>
  mutate(compromise = fct_recode(compromise,
                                 "Agree very strongly with 1" = "I agree strongly with A",
                                 "Agree with 1" = "I agree with A",
                                 "Agree with 2" = "I agree with B",
                                 "Agree very strongly with 2" = "I agree strongly with B",
                                 "Agree with neither" = "I agree with neither")) |>
  mutate(diverse = factor(diverse, levels = c("I agree strongly with A",
                                              "I agree with A",
                                              "I agree with B",
                                              "I agree strongly with B",
                                              "I agree with neither"))) |>
  mutate(diverse = fct_recode(diverse,
                              "Agree very strongly with 1" = "I agree strongly with A",
                              "Agree with 1" = "I agree with A",
                              "Agree with 2" = "I agree with B",
                              "Agree very strongly with 2" = "I agree strongly with B",
                              "Agree with neither" = "I agree with neither")) |>
  mutate(federalism = factor(federalism, levels = c("I agree strongly with A",
                                                    "I agree with A",
                                                    "I agree with B",
                                                    "I agree strongly with B",
                                                    "I agree with neither"))) |>
  mutate(federalism = fct_recode(federalism,
                                 "Agree very strongly with 1" = "I agree strongly with A",
                                 "Agree with 1" = "I agree with A",
                                 "Agree with 2" = "I agree with B",
                                 "Agree very strongly with 2" = "I agree strongly with B",
                                 "Agree with neither" = "I agree with neither")) |>
  mutate(attend_comm = as.double(attend_comm)) |>
  mutate(attend_protest = as.double(attend_protest))

# Combine
data_comp <- bind_rows(afro_comp, baseline_comp) |>
  mutate(source = ifelse(row_number() <= nrow(afro_comp), "Afrobarometer", "Study")) |>
  mutate(source = factor(source, levels = c("Afrobarometer", "Study"))) |>
  mutate(withinwt_hh = ifelse(source == "Study", 1, withinwt_hh))

# Participation comparison
weighted_averages <- data_comp %>%
  select(source, withinwt_hh, attend_protest, attend_comm) %>%
  pivot_longer(cols = c(attend_protest, attend_comm),
               names_to = "activity", values_to = "value") %>%
  group_by(source, activity) %>%
  summarize(
    weighted_avg = weighted.mean(value, w = withinwt_hh, na.rm = TRUE),
    .groups = "drop"
  )

weighted_averages %>%
  ggplot(aes(x = weighted_avg, y = activity, fill = source)) +
  geom_col(position = "dodge") +
  scale_x_continuous(labels = percent_format(accuracy = 1)) +
  scale_y_discrete(
    labels = c("attend_protest" = "Protest Attendance",
               "attend_comm" = "Community Meeting Attendance")) +
  labs(x = "Weighted Participation Rate", y = NULL,
       fill = "Source:", caption = "Source: Afrobarometer (2020)") +
  theme_minimal() +
  theme(legend.position = "right", strip.text = element_text(hjust = 0.5))

ggsave_both(name = "participation_ab_baseline_compari", width = 10, height = 6)
```

### Figure A3: Participation by Group

Source: `07_afrobarometer group comparison.R`. Participation index across demographic groups.

```{r figA3-participation-groups}
#| label: fig-participation-groups
#| fig-width: 10
#| fig-height: 8

participation_summary <- clean_ab |>
  summarise(
    pct_college = weighted.mean(participation[college == 1], w = withinwt_hh[college == 1], na.rm = TRUE),
    pct_young_25 = weighted.mean(participation[young_25 == 1], w = withinwt_hh[young_25 == 1], na.rm = TRUE),
    pct_all = weighted.mean(participation, w = withinwt_hh, na.rm = TRUE)
  ) |>
  pivot_longer(cols = everything(), names_to = "group", values_to = "pct") |>
  mutate(
    group = factor(group, levels = c("pct_college", "pct_young_25", "pct_all"),
                   labels = c("College (Graudate or Enrolled)", "Young (Under 25)", "Whole Sample"))
  )

ggplot(participation_summary, aes(x = pct, y = group)) +
  geom_col() +
  labs(title = "Participation by Group", y = NULL,
       x = "Average Participation Index\n (Community Attendence, Protest, Raise Issue)",
       subtitle = "Source: Afrobarometer, Ethiopia (2020).",
       caption = "Note: Means use sampling weights.") +
  scale_y_discrete(
    labels = c("College (Graudate or Enrolled)" = "College (Graduate or Enrolled)\n n = 127",
               "Young (Under 25)" = "Young (Under 25)\n n = 640",
               "Whole Sample" = "Whole Sample\n n = 2378")) +
  theme_minimal()

ggsave_both(name = "participation_groups", width = 10, height = 8)
```

### Figure A4: Protest by Group

Source: `07_afrobarometer group comparison.R`. Protest participation rates across demographic groups.

```{r figA4-protest-groups}
#| label: fig-protest-groups
#| fig-width: 10
#| fig-height: 8

protest_summary <- clean_ab |>
  summarise(
    pct_college = weighted.mean(protest[college == 1], w = withinwt_hh[college == 1], na.rm = TRUE),
    pct_young_25 = weighted.mean(protest[young_25 == 1], w = withinwt_hh[young_25 == 1], na.rm = TRUE),
    pct_all = weighted.mean(protest, w = withinwt_hh, na.rm = TRUE)
  ) |>
  pivot_longer(cols = everything(), names_to = "group", values_to = "pct") |>
  mutate(
    group = factor(group, levels = c("pct_college", "pct_young_25", "pct_all"),
                   labels = c("College (Graudate or Enrolled)", "Young (Under 25)", "Whole Sample"))
  )

ggplot(protest_summary, aes(x = pct, y = group)) +
  geom_col() +
  labs(title = "Percent Participated in Protest by Group", y = NULL,
       x = "% Protested",
       subtitle = "Source: Afrobarometer, Ethiopia (2020).",
       caption = "Note: Protest rates use sampling weights.") +
  scale_y_discrete(
    labels = c("College (Graudate or Enrolled)" = "College (Graduate or Enrolled)\n n = 127",
               "Young (Under 25)" = "Young (Under 25)\n n = 640",
               "Whole Sample" = "Whole Sample\n n = 2378")) +
  scale_x_continuous(labels = percent_format(accuracy = 1)) +
  theme_minimal()

ggsave_both(name = "protest_groups", width = 10, height = 8)
```

### Figure A5: Voting by Group

Source: `07_afrobarometer group comparison.R`. Voting rates across demographic groups.

```{r figA5-voting-groups}
#| label: fig-voting-groups
#| fig-width: 10
#| fig-height: 8

voting_summary <- clean_ab |>
  summarise(
    pct_college = weighted.mean(vote[college == 1], w = withinwt_hh[college == 1], na.rm = TRUE),
    pct_young_25 = weighted.mean(vote[young_25 == 1], w = withinwt_hh[young_25 == 1], na.rm = TRUE),
    pct_all = weighted.mean(vote, w = withinwt_hh, na.rm = TRUE)
  ) |>
  pivot_longer(cols = everything(), names_to = "group", values_to = "pct") |>
  mutate(
    group = factor(group, levels = c("pct_college", "pct_young_25", "pct_all"),
                   labels = c("College (Graudate or Enrolled)", "Young (Under 25)", "Whole Sample"))
  )

ggplot(voting_summary, aes(x = pct, y = group)) +
  geom_col() +
  labs(title = "Percent Voted by Group", y = "Group", x = "% Voted",
       subtitle = "Source: Afrobarometer, Ethiopia (2020).",
       caption = "Note: Voting rates use sampling weights.") +
  scale_y_discrete(
    labels = c("College (Graudate or Enrolled)" = "College (Graduate or Enrolled)\n n = 127",
               "Young (Under 25)" = "Young (Under 25)\n n = 640",
               "Whole Sample" = "Whole Sample\n n = 2378")) +
  scale_x_continuous(labels = percent_format(accuracy = 1)) +
  theme_minimal()

ggsave_both(name = "voting_groups", width = 10, height = 8)
```

### Figure A6: Political Views Comparison

Source: `08_compare_afro_baseline.R`. Comparison of political attitudes between Afrobarometer and study sample.

```{r figA6-views}
#| label: fig-views
#| fig-width: 15
#| fig-height: 6

cross_tabs <- data_comp |>
  select(source, withinwt_hh, compromise, diverse, federalism) |>
  pivot_longer(-c(source, withinwt_hh), names_to = "views", values_to = "value") |>
  mutate(value = fct_collapse(value,
                              "Agree with 1" = c("Agree very strongly with 1", "Agree with 1"),
                              "Neither" = c("Agree with neither"),
                              "Agree with 2" = c("Agree with 2", "Agree very strongly with 2")),
         value = factor(value,
                        levels = c("Agree with 1", "Neither", "Agree with 2"),
                        ordered = TRUE),
         value = fct_rev(value)) |>
  group_by(source, views) |>
  count(value, wt = withinwt_hh) |>
  filter(!value %in% c("Don't know", "Refused")) |>
  drop_na()

cross_tabs <- cross_tabs |>
  mutate(views = case_when(
    views == "compromise" ~ "1: No political compromise\n2: Political compromise",
    views == "diverse" ~ "1: Diversity is strength\n2: Homogeneity is strength",
    views == "trust" ~ "1: Most can be trusted\n2: Be wary of others",
    views == "federalism" ~ "1: Keep ethnic federalism\n2:change to geographic federalism")) |>
  group_by(source, views) |>
  mutate(percent = n / sum(n) * 100)

cross_tabs |>
  ggplot(aes(y = value, x = percent, fill = source)) +
  geom_col(position = "dodge") +
  facet_wrap(vars(views)) +
  scale_x_continuous(labels = scales::percent_format(scale = 1)) +
  labs(x = "Percent of respondents",
       y = "views", fill = "Compare:",
       caption = "Source: Afrobarometer (2020)")

ggsave_both(name = "views", width = 15, height = 6)
```

## Other Results

### Figure A7: Conjoint by Public Sector Plans

Source: `06-decentralize-cjt.R`. Marginal means split by post-graduation public sector plans.

```{r figA7-public}
#| label: fig-cjt-public
#| fig-width: 8
#| fig-height: 7

fed_cjt <- fed_cjt |>
  mutate(public_plans = ifelse(q26 == "Public sector or civil service (ex. public university, government ministry)", "Public sector", "Other")) |>
  mutate(public_plans = factor(public_plans, levels = c("Other", "Public sector")))

mm_public <- fed_cjt %>%
  select(outcome, hiring, police, culture, borders, public_plans) %>%
  drop_na() %>%
  cj(data = .,
     formula = outcome ~ hiring + police + culture + borders,
     by = ~ public_plans,
     estimate = "mm")

mm_public |>
  left_join(features) %>%
  left_join(levels) %>%
  ggplot(aes(x = estimate, y = clean_level,
             xmin = lower, xmax = upper,
             color = public_plans)) +
  geom_pointrange(fatten = 3, linewidth = 2, shape = 21, fill = "white",
                  position = position_dodge(width = 0.5)) +
  theme(legend.position = "top", strip.text = element_text(hjust = 0.5)) +
  labs(x = "Marginal Means\nPr(Profile chosen | attribute)",
       y = NULL,
       color = "Post-graduation plans:") +
  facet_wrap(vars(clean_feature), ncol = 1, scales = "free_y") +
  scale_x_continuous(labels = scales::percent) +
  scale_y_discrete(labels = scales::label_wrap(20)) +
  scale_color_manual(values = c(blue, red)) +
  geom_vline(xintercept = 0.5, lty = 2)

ggsave_both(name = "fed_cjt_public", width = 8)
```

### Table A1: Post-Hoc Power Analysis

Source: `05-vignette-experiment.R`. Minimum detectable effect sizes for pairwise comparisons.

```{r tabA1-posthoc}
#| label: tab-posthoc

n_control <- baseline |> filter(federalism_treat == "Control (no info)") |> nrow()
n_t1 <- baseline |> filter(federalism_treat == "R is in majority") |> nrow()
n_t2 <- baseline |> filter(federalism_treat == "R is in minority") |> nrow()
n_t3 <- baseline |> filter(federalism_treat == "R is in majority, but most co-ethnics in minority") |> nrow()
alpha <- 0.05
power <- 0.80

min_detectable_d <- function(n1, n2, alpha = 0.05, power = 0.80) {
  pwr.t2n.test(n1 = n1, n2 = n2, sig.level = alpha, power = power)$d
}

d_t1 <- min_detectable_d(n_control, n_t1, alpha, power)
d_t2 <- min_detectable_d(n_control, n_t2, alpha, power)
d_t3 <- min_detectable_d(n_control, n_t3, alpha, power)

tab <- tribble(~Test, ~`Min. Detectable d`,
               "Control vs. T1", d_t1,
               "Control vs. T2", d_t2,
               "Control vs. T3", d_t3) |>
  mutate(across(c(`Min. Detectable d`), ~ round(., 3)))

knitr::kable(tab, caption = "Table A1: Minimum detectable Cohen's d for pairwise comparisons.")

# Also export LaTeX
invisible(capture.output(
  stargazer::stargazer(tab, summary = FALSE,
                       title = "Minimum detectable Cohen's d for pairwise comparisons.",
                       label = "post-hoc",
                       out = "figures/post-hoc.tex",
                       style = "ajps",
                       rownames = FALSE)
))
```

### Table A2: Conjoint Marginal Means

Source: `06-decentralize-cjt.R`. Marginal means with p-values for conjoint attributes.

```{r tabA2-cjt-mms}
#| label: tab-cjt-mms

table_data <- mms %>%
  select(feature, level, estimate, p_value) %>%
  rename(Feature = feature, Level = level,
         Coefficient = estimate, `P-value` = p_value) %>%
  arrange(Feature, Level)

table_data$Level <- c("Central Gov Hires", "States Hire", "State Police", "National Police",
                       "State Language and Culture", "Central Gov Language and Culture",
                       "Geography", "Ethnic Group")

knitr::kable(table_data, digits = 3, caption = "Table A2: Conjoint Analysis Results (MMs)")

invisible(capture.output(
  stargazer(as.data.frame(table_data), summary = FALSE,
            title = "Conjoint Analysis Results (MMs)",
            rownames = FALSE, digits = 3,
            label = "tab:federal_cjt_mms",
            out = "figures/federal_cjt_mms.tex")
))
```

### Figure A8: Conjoint AMCEs by Ethnicity

Source: `06-decentralize-cjt.R`. AMCEs split by three ethnic groups (Oromo, Amhara, Other).

```{r figA8-ethnic3-amce}
#| label: fig-cjt-ethnic3-amce
#| fig-width: 8
#| fig-height: 7

mm_ethnicity3_amce <- fed_cjt %>%
  mutate(three_ethn = as.factor(three_ethn)) |>
  select(outcome, hiring, police, culture, borders, three_ethn) %>%
  drop_na() %>%
  cj(data = .,
     formula = outcome ~ hiring + police + culture + borders,
     by = ~ three_ethn,
     estimate = "amce")

mm_ethnicity3_amce |>
  left_join(features) %>%
  left_join(levels) %>%
  ggplot(aes(x = estimate, y = clean_level,
             xmin = lower, xmax = upper,
             color = three_ethn)) +
  geom_pointrange(fatten = 3, linewidth = 2, shape = 21, fill = "white",
                  position = position_dodge(width = 0.5)) +
  theme(legend.position = "top", strip.text = element_text(hjust = 0.5)) +
  labs(x = "Average marginal component effects (AMCEs)\nPr(Profile chosen | attribute) - Pr(Profile chosen | baseline)",
       y = NULL,
       color = "Ethnic Group",
       caption = "Baseline attribute centered at zero.") +
  facet_wrap(vars(clean_feature), ncol = 1, scales = "free_y") +
  scale_x_continuous(labels = scales::percent) +
  scale_y_discrete(labels = scales::label_wrap(20)) +
  geom_vline(xintercept = 0, lty = 2) +
  scale_color_manual(values = c(red, green, blue))

ggsave_both(name = "fed_cjt_ethnic3_amce", width = 8)
```

### Table A3: Main Results

Source: `05-vignette-experiment.R`. OLS regression of federalism support on treatment conditions.

```{r tabA3-main}
#| label: tab-main
#| results: asis

mod <- lm(q45_1 ~ federalism_treat, data = baseline)

stargazer::stargazer(mod,
                     title = "Table A3: Federalism Vignettes Main Results",
                     dep.var.labels = "Support for Federal Power",
                     style = "ajps",
                     type = "html",
                     covariate.labels = c("Majority", "Majority, Coethnics Minority", "Minority"),
                     label = "tab:main-result",
                     out = "figures/main-table.tex")
```

### Table A4: Co-Ethnic President Interaction

Source: `05-vignette-experiment.R`. OLS regression interacting treatment with co-ethnic president indicator.

```{r tabA4-pres}
#| label: tab-pres
#| results: asis

mod_pres <- lm(q45_1 ~ federalism_treat * is_oromo, data = baseline)

stargazer::stargazer(mod_pres,
                     title = "Table A4: Federalism Vignettes Co-Ethnic President Results",
                     dep.var.labels = "Support for Federal Power",
                     style = "ajps",
                     type = "html",
                     covariate.labels = c("Majority", "Minority", "Majority, Coethnics Minority", "Oromo",
                                          "Majority*Oromo", "Minority*Oromo", "Majority, Coethnics Minority*Oromo"),
                     label = "tab:pres-result",
                     out = "figures/pres-table.tex")
```

### Table A5: One-Party Rule Interaction

Source: `05-vignette-experiment.R`. OLS regression interacting treatment with one-party rule support.

```{r tabA5-1party}
#| label: tab-1party
#| results: asis

mod_1p <- lm(q45_1 ~ federalism_treat * one_party, data = baseline)

stargazer::stargazer(mod_1p,
                     title = "Table A5: Federalism Vignettes One-Party Rule Results",
                     dep.var.labels = "Support for Federal Power",
                     style = "ajps",
                     type = "html",
                     covariate.labels = c("Majority", "Majority, Coethnics Minority", "Minority",
                                          "Supports 1-Party Rule", "Majority*1-Party",
                                          "Majority, Coethnics Minority*1-Party", "Minority*1-Party"),
                     label = "tab:1-party-result",
                     out = "figures/1-party-table.tex")
```

### Table A6: Tolerance Interaction

Source: `05-vignette-experiment.R`. OLS regression interacting treatment with ethnic tolerance.

```{r tabA6-tolerance}
#| label: tab-tolerance
#| results: asis

mod_tol <- lm(q45_1 ~ federalism_treat * tolerance, data = baseline)

stargazer::stargazer(mod_tol,
                     title = "Table A6: Federalism Vignettes Tolerance Results",
                     dep.var.labels = "Support for Federal Power",
                     style = "ajps",
                     type = "html",
                     covariate.labels = c("Majority", "Majority, Coethnics Minority", "Minority",
                                          "Ethnic Tolerance", "Majority*Tolerance",
                                          "Majority, Coethnics Minority*Tolerance", "Minority*Tolerance"),
                     label = "tab:tolerance-result",
                     out = "figures/tolerance-table.tex")
```

### Figure A9: Conjoint AMCEs (Overall)

Source: `06-decentralize-cjt.R`. Overall average marginal component effects.

```{r figA9-amces}
#| label: fig-cjt-amces
#| fig-width: 8
#| fig-height: 7

amces <- cj(data = fed_cjt,
            formula = form,
            id = ~ r_id,
            estimate = "amce") |>
  as_tibble()

amces %>%
  left_join(features) %>%
  left_join(levels) %>%
  ggplot(
    aes(x = estimate, y = clean_level,
        xmin = lower, xmax = upper)) +
  geom_pointrange(fatten = 3, linewidth = 2, shape = 21, fill = "white", color = blue) +
  theme(legend.position = "none", strip.text = element_text(hjust = 0.5)) +
  labs(x = "Average marginal component effects (AMCEs)\nPr(Profile chosen | attribute) - Pr(Profile chosen | baseline)",
       y = NULL) +
  facet_wrap(vars(clean_feature), ncol = 1, scales = "free_y") +
  scale_x_continuous(labels = scales::percent) +
  scale_y_discrete(labels = scales::label_wrap(30)) +
  scale_color_manual(values = palette[c(2, 4, 6, 8)]) +
  geom_vline(xintercept = 0, lty = 2)

ggsave_both(name = "federal_cjt_amces", width = 8)
```

### Figure A10: AMCE Differences by Border Type

Source: `06-decentralize-cjt.R`. AMCEs conditional on border type (geographic vs. ethnic) and their difference.

```{r figA10-amce-diff-borders}
#| label: fig-amce-diff-borders
#| fig-width: 10
#| fig-height: 7

amces_borders <- fed_cjt |>
  select(outcome, hiring, police, culture, borders) |>
  drop_na() |>
  cj(formula = outcome ~ hiring + police + culture,
     by = ~ borders,
     estimate = "amce") |>
  left_join(features) |>
  left_join(levels) |>
  as_tibble()

amces_borders_diff <- fed_cjt |>
  select(outcome, hiring, police, culture, borders) |>
  drop_na() |>
  cj(formula = outcome ~ hiring + police + culture,
     by = ~ borders,
     estimate = "amce_diff") |>
  left_join(features) |>
  left_join(levels) |>
  as_tibble()

rbind(amces_borders, amces_borders_diff) |>
  mutate(clean_feature = case_when(str_detect(clean_feature, "police") ~ "Policing",
                                   str_detect(clean_feature, "language") ~ "Cultural policy",
                                   str_detect(clean_feature, "hiring") ~ "Public hiring")) |>
  mutate(BY = case_when(BY == "State borders are drawn according to geography" ~ "Geographic borders",
                        BY == "State borders are drawn according to majority ethnic group in the area" ~ "Ethnic borders",
                        .default = "Difference in AMCEs"),
         BY = factor(BY, levels = c("Geographic borders", "Ethnic borders", "Difference in AMCEs"))) |>
  ggplot(aes(x = estimate, y = clean_level,
             xmin = lower, xmax = upper)) +
  geom_pointrange(fatten = 3, linewidth = 2, shape = 21, fill = "white", color = blue) +
  theme(legend.position = "none", strip.text = element_text(hjust = 0.5)) +
  labs(x = "Average marginal component effects (AMCEs)\nPr(Profile chosen | attribute) - Pr(Profile chosen | baseline)",
       y = NULL) +
  facet_grid(rows = vars(clean_feature),
             cols = vars(BY), scales = "free_y") +
  scale_x_continuous(labels = scales::percent) +
  scale_y_discrete(labels = scales::label_wrap(30)) +
  geom_vline(xintercept = 0, lty = 2)

ggsave_both(name = "amce_diff_borders", width = 10)
```

### Table A7: Ethnic Violence Difference in Marginal Means

Source: `06-decentralize-cjt.R`. Difference in marginal means between ethnic violence supporters and opponents.

```{r tabA7-eth-violence-dmms}
#| label: tab-eth-violence-dmms

dmm_eth_violence <- fed_cjt %>%
  select(outcome, hiring, police, culture, borders, ev_dichotomous) %>%
  drop_na() %>%
  cj(data = .,
     formula = outcome ~ hiring + police + culture + borders,
     by = ~ ev_dichotomous,
     estimate = "mm_differences") %>%
  mutate(ev_dichotomous = factor(ev_dichotomous, levels = c(0, 1),
                                 labels = c("No", "Yes"))) %>%
  mutate(feature = stringi::stri_trans_totitle(feature),
         feature = factor(feature, levels = c("Hiring", "Police", "Culture", "Borders")))

table_data_ev <- dmm_eth_violence %>%
  select(feature, level, estimate, p) %>%
  rename(Feature = feature, Level = level, DMM = estimate) %>%
  arrange(Feature)

table_data_ev$Level <- c("Central Gov Hires", "States Hire", "State Police", "National Police",
                          "State Language and Culture", "Central Gov Language and Culture",
                          "Geography", "Ethnic Group")

knitr::kable(table_data_ev, digits = 3,
             caption = "Table A7: Difference in Marginal Means (comparing supporters of ethnic violence against opponents)")

invisible(capture.output(
  stargazer(as.data.frame(table_data_ev), summary = FALSE,
            title = "Difference in Marginal Means (comparing supporters of ethnic violence against opponents)",
            rownames = FALSE, digits = 3,
            label = "tab:eth_violence_cjt_dmms",
            out = "figures/eth_violence_cjt_dmms.tex")
))
```

### Table A8: Marginal Means by Ethnicity (3 Groups)

Source: `06-decentralize-cjt.R`. Marginal means for each conjoint condition by ethnic group.

```{r tabA8-3-ethnicities}
#| label: tab-3-ethnicities

mm_ethnicity3 <- fed_cjt %>%
  mutate(three_ethn = as.factor(three_ethn)) |>
  select(outcome, hiring, police, culture, borders, three_ethn) %>%
  drop_na() %>%
  cj(data = .,
     formula = outcome ~ hiring + police + culture + borders,
     by = ~ three_ethn,
     estimate = "mm")

mm_ethnicity3$condition <- c("Central Gov", "State", "State", "Central Gov",
                              "State", "Central Gov", "Geography", "Ethnicity",
                              "Central Gov", "State", "State", "Central Gov",
                              "State", "Central Gov", "Geography", "Ethnicity",
                              "Central Gov", "State", "State", "Central Gov",
                              "State", "Central Gov", "Geography", "Ethnicity")

table_data_eth <- mm_ethnicity3 %>%
  select(feature, condition, three_ethn, estimate, lower, upper) %>%
  arrange(feature, three_ethn) %>%
  mutate(estimate = round(estimate, 3),
         lower = round(lower, 3),
         upper = round(upper, 3)) %>%
  rename(Feature = feature, `Ethnic Group` = three_ethn,
         `Marginal Mean` = estimate, Lower = lower, Upper = upper)

knitr::kable(table_data_eth, caption = "Table A8: Marginal Means for each condition based on Ethnic Group")

invisible(capture.output(
  stargazer(table_data_eth, title = "Marginal Means for each condition based on Ethnic Group",
            summary = FALSE, rownames = FALSE,
            label = "tab:3_ethnicities_cjt_mms",
            out = "figures/3_ethnicities_cjt_mms.tex")
))
```

## Balance Tests

### Figure A11: Conjoint Attribute Frequencies

Source: `06-decentralize-cjt.R`. Diagnostic check that conjoint attributes are uniformly distributed.

```{r figA11-frequencies}
#| label: fig-cjt-frequencies
#| fig-width: 8
#| fig-height: 7

plot(cj_freqs(fed_cjt, form, id = ~r_id))

ggsave_both(name = "fed_cjt_frequencies", width = 8)
```

### Figure A12: Conjoint Balance on Multi-Ethnic Status

Source: `06-decentralize-cjt.R`. Balance test showing conjoint attributes are uncorrelated with multi-ethnic identity.

```{r figA12-balance}
#| label: fig-cjt-balance
#| fig-width: 8
#| fig-height: 7

fed_cjt <- fed_cjt %>%
  mutate(multiethnic_bin = case_when(multiethnic == TRUE ~ 1,
                                     multiethnic == FALSE ~ 0))
plot(mm(fed_cjt, multiethnic_bin ~ hiring + police + culture + borders, id = ~r_id),
     vline = mean(fed_cjt$multiethnic_bin, na.rm = TRUE))

ggsave_both(name = "fed_cjt_multiethnic_balance", width = 8)
```

### Figure A13: Conjoint Carryover Test

Source: `06-decentralize-cjt.R`. Tests for carryover effects between conjoint profiles.

```{r figA13-carryover}
#| label: fig-cjt-carryover
#| fig-width: 8
#| fig-height: 7

fed_cjt$profile <- factor(fed_cjt$profile)
plot(cj(fed_cjt, form, id = ~r_id, by = ~profile, estimate = "mm"), group = "profile", vline = 0.5)

ggsave_both(name = "fed_cjt_carry_over", width = 8)
```
